Comparison of Multibody Dynamics Solver Performance: Synthetic versus Realistic Data

In the area of robotics simulation, multibody dynamics plays an important role in designing and controlling robots, especially when the robot contacts the environment. Contacts give rise to non-penetration and friction constraints, which are nonsmooth and nonlinear. One way to simulate such systems is through the use of a discrete-time multibody dynamics model in the form of a nonlinear complementarity problem (NCP), for which, finding a solution is known to be NP-hard [1]. In situations where analytical solutions don’t exist, a suite of numerical solutions accessible through a benchmarking framework is useful to fairly evaluate performance of different computer algorithms. However, many algorithm designers don’t have easy access to test data from physical simulators. Under such circumstances, randomized data are used to test the performance of solution algorithms. In this paper, we present our Benchmark Problems of Multibody Dynamics (BPMD) framework and database, with the data sets from different physics engines, as a benchmarking platform. Then we compare the performance of several solvers on synthetic data and simulation data, to show the superiority of testing solution algorithms with simulation data over testing with synthetic data. We will show that algorithm tested only on synthetic data often fail to solve problem obtained from physics simulations, to demonstrate the benefit of BPMD database. ∗Address all correspondence to this author. INTRODUCTION The instantaneous-time equations of motion for multibody dynamic systems with unilateral contacts take the form of a differential complementarity problem (dCP). This equation must be integrated to simulate the system, but this task is made difficult by the fact that the dCP is nonlinear and nonsmooth.1 The class of integration methods known as “time-stepping methods,” first developed by Moreau, discretize the dCP with respect to time, converting the integration problem into one of solving a series of complementarity problems (CPs) [2–4]. The robustness and efficiency of CP solution algorithms are significant not only to multibody simulation, but also to modelpredictive robot control when the robot performs tasks involving intermittent contact. In both simulation and control, solving the CPs quickly and reliably is important. Unfortunately, the CPs have poor mathematical properties, including singularities, nonuniqueness, and NP-hardness [1]. These properties make it difficult to design solution algorithms guaranteed to converge to solutions of the CPs. However, a linearized version of the problem leads to a series of linear complementarity problems (LCPs) with good properties: a solution exists and can be found by Lemke’s algorithm.2 Nonetheless, there is still no agreed upon “best” algorithm for solving the CPs that arise in multibody dynamics 1Nonsmoothness is an artifact of imposing non-penetration between rigid bodies and modelling stick-slip friction precisely 2Admittedly, while the properties are better, the problem size is larger, which can increase solution time [5] 1 Copyright c © 2015 by ASME problems. These facts have motivated researchers to develop many different algorithms, primarily taking one of the following two approaches: 1. Modify the model, via linearization, relaxation, or regularization, to improve the mathematical properties and make them easier to solve with existing algorithms [6–8]. 2. Develop new algorithms to solve the model without modification. These fall into three main categories: pivoting methods, matrix-splitting methods, and second-order generalized Newton methods [9–11]. Whenever a multibody dynamics solver is proposed, it is desirable to test its performance and compare it with existing solvers. This should be done using data from actual multibody systems problems, such as robotic manipulation, mechanical assembly, and driving on realistic terrain, but most researchers do not have easy access to such data. Therefore the matrices of the test problems are constructed by a combination of pseudorandom number generation and matrix operations, so that the vectors and matrices of the problem satisfy certain properties known to be exhibited by physically-based problems. Such “synthetic” data may be useful in initial testing of solver accuracy and convergence rates [12, 13], but the results presented here show that solvers can behave quite differently (typically more poorly) when faced with physically-based problems. In this paper, we recommend that solvers be tested with physically-based data, and moreover, with data from a set of benchmark problems, so that all solvers can be compared fairly. We support this recommendation with a solver performance study that demonstrates significant differences in solver performance when applied to synthetic and physically-based data. This study was made simple to conduct by the Benchmark Problems in Multibody Dynamics (BPMD) database [14], whose primary goal is to facilitate the fair and thorough comparison of solvers used in multibody dynamics simulation software. For the remainder of this paper, we limit our discussion to LCPs, where we construct the complementarity conditions for unilateral constraints to model contacts, and then linearize the friction cones to multi-faceted polyhedra. The structure of the paper follows: We first introduce the general form of a CP, a specific LCP that is commonly used in multibody dynamics software, and several popular solvers. Before discussing the results of the comparative study, the Benchmark Problems for Multibody Dynamics (BPMD) database and framework are introduced. Next we present the methods that are used in this paper: pivoting, matrix-splitting and second-order generalized Newton’s methods. These methods are applied to the simulation data, and the synthetic data with the same matrix size as those from simulation, to test their convergence rate features. The last section presents the comparison results of solver performance between the cases of synthetic and physically-based data. COMPLEMENTARITY PROBLEMS AND SOLVERS Complementarity Problem Given a known vector function f(x) ∈ <, find x ∈ < (where n is a positive integer) satisfying the three conditions: x ≥ 0, f(x) ≥ 0, and x f(x) = 0. For brevity, we use the following short-hand notation: 0 ≤ x ⊥ f(x) ≥ 0. (1) When f(x) is a linear function of x, then the problem is LCP. A rigid body dynamics model with contact interactions is comprised of the Newton-Euler equation that determine dynamic motions of the objects and a contact model that enforces nonpenetration between bodies and models physically correct dry friction behaviour. Specifically, the contact model must allow stick-slip transitions, ensure maximal energy dissipation at sliding contacts, and ensure that the contact impulses lie within their respective friction cones. With fixed time step, a linearized, discrete-time version of the Newton-Euler equation and contact model can be expressed as follows [3, 4]: Mν = Gnp `+1 n +Gfp `+1 f +Mν ` + pext (2) 00 0  ≤ ρ`+1 n ρ f σ  ⊥ p`+1 n p f s  ≥ 00